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ISOTHERMAL ELASTOHYDRODYNAMIC LUBRICATION OF POINT CONTACTS 

I - THEORETICAL FORMULATION 
by Bernard J. Hamrock and Duncan Dowson* 

Lewis Research Center 

SUMMARY 

The isothermal elastohydrodynamic lubrication (EHL) of a point contact was ana- 
lyzed numerically by simultaneously solving the elasticity and Reynolds equations. In 
the elasticity analysis the contact zone was divided into equal rectangular areas, and it 
was assumed that a uniform pressure was applied over each area. In the numerical 
analysis of the Reynolds equation, a phi analysis (where phi is equal to the pressure 
times the film thickness to the 3/2 power) was used to help the relaxation process. The 
EHL point contact analysis is applicable for the entire range of elliptical parameters and 
is valid for any combination of rolling and sliding within the contact. 


INTRODUCTION 

In many contacts between machine elements, forces are transmitted through thin, 
but' continuous, fluid films. These fluid films as related to hydrodynamic lubrication in 
journal and thrust bearings have been well understood for some time, and experimental 
work confirms the theory. In the early 1900's, it was recognized by Martin (ref. 1) that 
many loaded contacts of low geometrical conformity such as gears and rolling-contact 
bearings probably behaved as though they were hydrodynamically lubricated. As opposed 
to the journal and thrust bearing counterpart, the original hydrodynamic lubrication 
theory of gears and rolling -contact bearings differed substantially from experimental 
findings. Only in recent years has the consideration of elastic deformation of contacts 
been coupled to hydrodynamics to yield closer agreement of theory with experiments. 

Elastohydrodynamic lubrication (EHL) deals with the lubrication of elastic contacts. 
The analysis requires the simultaneous solution of the elasticity and Reynolds equations. 
The EHL theory differs from conventional hydrodynamic theory in the following ways: 

* Professor of Mechanical Engineering, Leeds University, Leeds, England. 



(1) In defining the film thickness in EHL theory, elastic deformation of the contact 
is considered. 

(2) The viscosity is no longer independent of pressure, as is frequently assumed in 
conventional hydrodynamic theory. 

Figure 1 shows the differences between the Martin (ref. 1), Hertzian, and EHL lub- 
rication and contact conditions. Pressure and film thickness curves are shown for the 
respective conditions, along with a formula for minimum film thickness. The preceding 
definition of EHL covers not only such highly stressed lubricated machine elements as 
ball and roller bearings and gears, but also the performance of soft rubber seals, bear- 
ings with soft liners, and human joints. One may, therefore, classify these two types 
of EHL as hard EHL (ball and roller bearings and gears) and soft EHL (rubber seals, 
bearings with soft liners, and human joints). In this report we are concerned only with 
hard EHL. 

When two solids are in contact under zero-load condition, one of two types of con- 
tact is experienced: 

(1) Point contact, in which the two solids touch at a single point, as in ball bearings 

(2) Line contact, in which the two solids touch along a straight or curved line, as in 

roller bearings 

After a load is applied, the point expands to an ellipse and the line to a rectangle. Al- 
though we are concerned with loaded contacts, it is convenient to distinguish between 
these situations by referring to them as being either point or line contact. 

Most of the work done to date in EHL has dealt with line contacts. Grubin (ref. 2) 
was the first to attempt a solution to the isothermal EHL line-contact problem. In Gru- 
bin's analysis it was assumed that the shape of the elastically deformed solids in a 
highly loaded lubricated contact was the same as the shape produced in a dry (Hertzian) 
contact. This assumption facilitated the solution of the Reynolds equation in the inlet 
region of the contact and enabled the separation of the solids in the central region of the 
contact to be determined with commendable accuracy. 

In references 3 and 4 an empirical formula for the isothermal line-contact EHL 
problem was obtained. This formula shows the effect of speed, load, and material prop- 
erties upon minimum film thickness and is based on their theoretical solutions. The ex- 
perimental observations concur with the minimum -film -thickness formula. Forty-five 
years transpired from Martin's (ref. 1) rigid isothermal hydrodynamic solution until 
Dowson and Higginson’s (refs. 3 and 4) minimum-film-thickness formula for an isother- 
mal EHL line contact was developed. 

The study of point contacts has mostly followed experimental lines. The first step 
toward a theoretical solution of the point- contact problem was presented by Archard and 
Cowking (ref. 5). They adopted an approach similar to that used by Grubin (ref. 2) for 
line-contact conditions. The Hertzian contact zone was assumed to form a parallel film 
region, and the generation of high pressure in the approaches to the Hertzian zone was 
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considered. Cheng (ref. 6) also used an approach similar to that used by Grubin (ref. 2) 
in determining a minimum film thickness for point contacts. 

The reason why a complete, isothermal, point-contact EHL solution has not yet 
been obtained is the extreme difficulty of the numerical coupling of the elasticity and 
Reynolds equations. Within the last year an interesting numerical solution to the point- 
contact EHL problem for a sphere near a plane has been presented by Ranger (ref. 7). 
This solution is presented in dimensional terms, which thereby limits its general usage. 
A puzzling feature of Ranger's work (ref. 7) is the fact that his resulting equation for 
the minimum film thickness has a positive load exponent, which contradicts experiments 
(e. g. , ref. 8). 

In the present report the theoretical foundation for the complete, isothermal, point- 
contact EHL problem is presented. The theory is applicable to varying values of ellip- 
tical parameters. Utilized in the theory is elasticity analysis developed by the authors 
in an earlier publication (ref. 9). In this analysis the contact zone was divided into 
equal rectangular areas, and a uniform pressure was applied over each area. The 
pressure-viscosity analysis of Roelands (ref. 10) is also used and is valid for any com- 
bination of rolling and sliding. Because of the length and complexity of the EHL theory, 
it is presented by itself. The results of the numerical analysis will appear in separate 
reports. 


GEOMETRY OF CONTACTING ELASTIC SOLIDS 

Two solids having different radii of curvature in a pair of principal planes (x and y) 
passing through the contact between the solids make contact at a single point under the 
condition of no applied load. Such a condition is called "point contact" and is shown in 
figure 2, where the radii of curvature are denoted by r's. In the analysis which follows, 
it is assumed that for convex surfaces (fig. 1) the curvature is positive but that for con- 
cave surfaces the curvature is negative. The curvature sum and difference are defined 
as 


where 



( 1 ) 

( 2 ) 
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_i_ _ i i 

R x r Ax r Bx 


( 3 ) 


1 _ 1 : 1 
R y r Ay r By 


(4) 


When the two solids in figure 2 have a normal load applied to them, the result is 
that the point expands to an ellipse, with a being the semimajor axis and b being the 
semiminor axis. The normal applied load F in figure 2 lies along the axis which 
passes through the center of the solids at the point of contact and is perpendicular to a 
plane which is tangential to both solids at the point of contact. For the special case 
where r^ = r^ and r B ^ = r Bx , the resulting contact is a circle rather than an 
ellipse. 

The elliptical eccentricity parameter k is defined as 


k 



(5) 


From reference 11, k can be written in terms of the curvature difference and the ellip- 
tical integrals of the first and second kind as 


where 


J(k) 


V 2jf- <?(i + r) 
<m - r) 



(6) 


(?) 


( 8 ) 


A one-point iteration method which has been used successfully in the past (ref. 12) is 
used, where 


k n+l = J ( k n> 


(9) 
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When k is known, the semimajor and semiminor axes of the contact ellipse can be 
written as 


I 


where 




and 

v Poisson ratio 

E modulus of elasticity 


( 10 ) 

(ID 

(12) 


REYNOLDS EQUATION 

The coordinate systems to be used are shown in figure 3. For the coordinate sys- 
tem (X, Y) the general Reynolds equations for EHL of point contacts from reference 13 
gives the following: 


+ -4f^ =-iph 

ax\i2rj ax/ dY\i2n aY/ ax 


8 „J V A + V B> 

+ — ^ ph 

3Y 2 


where 


u^ surface velocity of solid A in X direction 

r>j 

Ug surface velocity of solid B in X direction 

v^ surface velocity of solid A in Y direction 

Vg surface velocity of solid B in Y direction 

The surface velocities are assumed to be constant and, therefore, equation (13) can be 
rewritten as 
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(14) 


where 


-i/fihf i£) + 4 = 12 [u 

ax\ n ax/ 3Y\ v sy J [ ax 


■^=^+ V 


3Y 


u = 


U A + U B 


(15) 


v - 


V A + V B 


(16) 


By introducing V and 0, where 




2 2 
+ v 


0 = tan -1 /- 


© 


equation (14) becomes 


-4/fiii! 9 r\ +_9 s /fihf ajA = 12V r 

ax\ r? ax/ 3Y\ v ay/ [ 


cos 0 -^2^ + sin 0 -^2^ 
3X 3Y 


By letting 


x = * 


= JL 

b 

7 0 

/V 

y=y 

H 

_ h 

a 


R 

P=£- 

P 

= JL 

2q 


E’ 




equation (19) can be written as 


-i/fig! J?) + J. ±(e£ Jp) = i2uA\L s s i®. + *1*1 &. 

3 X\ r? 3X/ k 2 3Y\ r? 3Y/ \r/[ 3X k 3Y . 


(17) 

(18) 


(19) 


( 20 ) 


( 21 ) 
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where 


V V 0 

u = ? (22) 

RE’ 

Equation ( 21 ) is the Reynolds equation for which the dimensionless pressure P will 
be determined. However, before proceeding, the dimensionless density p, the dimen- 
sionless viscosity 77, and the dimensionless film thickness H need to be expressed. 


Density 

From reference 14 the dimensionless density for mineral oil can be written as 


p = i+ °- 009 p 

1 + 0. 026 p 


where 

9 

p gage pressure, ton/in. 

Therefore, the general expression for the dimensionless density can be written as 


P = 1 + 


«PE’ 

1 +/3PE' 


(23) 


where 


a,f 3 constants dependent on fluid 


Viscosity 

Roelands (ref. 10 ) has undertaken a wide-ranging study of the effect of pressure 
upon the viscosity of lubricants. For isothermal conditions, his equation (p. 95 , 
ref. 10) can be written as 


log 77 + I. 200 = (log 7j + 1. 200) 



( 24 ) 


where 
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p gage pressure, kgf/cm^ 

ri 0 atmospheric viscosity 

z viscosity pressure index, a dimensionless constant 

Rearranging terms gives 



(25) 


where 

Voo 0. 00631 N-sec/m 2 (0. 0631 cP) 
y constant equal to 19 609 N/cm 2 (28 440 lb/in. 2 ) 

In equations (23) and (25), care must be taken to ensure that the same dimensions are 
used in defining the constants. 


Film Thickness 

A simplified expression for defining the film thickness may be written as 


h(X, Y) = h Q + S(X, Y) + w(X, Y) 


(26) 


where 

h central film thickness 

o 

S(X, Y) separation due to geometry of ellipsoidal solids 
w(X, Y) elastic deformation 

The separation due to the geometry of the two ellipsoidal solids shown in figure 2 
can be adequately described by an equivalent ellipsoidal solid near a plane. The geomet- 
rical requirement is that the separation of the ellipsoidal solids in the initial and equiv- 
alent situations should be the same at equal values of X. Therefore, from figure 3 the 
separation due to the geometry of the two ellipsoids shown in figure 2 can be written for 
an ellipsoidal solid near a plane as 
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(27) 


S(X, Y) = 


O 

(X - mbr 

2R x 


, (Y - l a) 2 
2R„ 


where 

m constant used to determine length of inlet region 
l constant used to determine length of side leakage region 

For purposes of illustration we will continue to use the mesh described in figure 3, 
where m = 4 and 1=2. However, the equations developed will be written in general 
terms. 

Figure 4 gives a physical description of the film thickness and its components for an 
ellipsoidal solid near a plane. Here it is assumed that y = 0, or Y = la., in figure 3. 

Substituting equation (27) into equation (26) while at the same time making this equa- 
tion dimensionless gives 


H(X, Y) = H + + w(X, Y) 

0 2RR 2RR R 

x y 


(28) 


where 


R 

In reference 9 the elastic deformation (w in eq. (28)) of an equivalent ellipsoidal 
solid near a plane in contact and subjected to a Hertzian stress distribution has been 
evaluated numerically. In the analysis of reference 9 the contact zone was divided into 
equal rectangular areas, and it was assumed that a uniform pressure was applied over 
each area. The results of reference 9 indicate that the division of the semimajor and 
semiminor axes into five equal subdivisions is sufficient to obtain the elastic deforma- 
tion with adequate accuracy. Making use of the results of reference 9, figure 3 shows 
an example of dividing the area in and around the contact zone into equal rectangular 
areas. In figure 3 the outlet region of the contact is not extended as much as the inlet 
region since previous investigators (e. g. , ref. 15) have found the pressure gradients at 
the outlet to be very large, with the pressure quickly setting to zero. 

Therefore, by using figure 3 and the results of reference 9, the elastic deformation 
can be written as 


9 



( 29 ) 



n constant used to determine length of outlet region 
c number of equal divisions in semimajor axis 

d number of equal divisions in semiminor axis 

m = | k - i | +1 


n = 1 1 - j | + 1 



(In fig. 3, c = d = 5. ) Equation (31) points out more explicitly the meaning of equa- 
tion (29) while making use of figure 3. The elastic deformation at the center of the rec- 
tangular area Wg ^ (fig. 3) caused by the pressure on the various rectangular areas in 
and around the contact ellipse can be written as 
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P l, 1 °9, 5 + P 2, 1 °8, 5 + • • • + P 35? i D 2 7, 5 


Wq _ = 2 + P l,2 D 9,4 + P 2,2 D 8,4 + ’ * ' + P 35, 2 °27, 4 
9 ’ 5 n 


L + P l, 20 D 9, 16 + P 2, 20°8, 16 + 1 ■ * + P 35, 20 D 27, 16. 


PHI ( cp ) SOLUTION 

Having defined the density, viscosity, and film thickness we can return to the prob- 
lem of solving the Reynolds equation. It is well known (e. g. , ref. 16) that the dimen- 
sionless pressure P of the Reynolds equation plotted as a function of X exhibits a very 
localized pressure field, giving high values of 3P/3X and 3 2 P/3X 2 . Such a condition 
with high gradients is not welcomed when performing a numerical analysis by means of 
relaxation methods. In order to produce a much gentler curve, a parameter <p is in- 
troduced, where 

4> = PH 3/2 (32) 

The pressure P is small at large values of film thickness H and conversely. This 
substitution also has the advantage of eliminating all terms containing derivatives of 
products of H and P or H and <fi. Therefore, by using equation (32), equation (21) 
can be written as 

JLQL/h 3 / 2 - 3 0H 1 / 2 i»\l + J_ _3_f |/ h 3 / 2 - 3 ah 1 / 2 9H\ 

axLr?\ ax 2 ax/J k 2 aY|_77 \ ay 2 ay/_ 


= 12u/— \ [cos 9 ifr H > + sine 

\R/ L k ay _ 

Expanding this equation yields 
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Finite Difference Representation 

The finite dif ference method will be used to develop the various terms in equa- 
tion (33). Figure 5 shows the mesh to be used as related to the dimensionless coordi- 
nates X and Y. Equation (33) must be written for the point (i, j) in figure 5 by sub- 
stituting for the derivatives expressions that involve values of $, H, p, and rj at the 
surrounding points. 

At the three points X. * -, X, ., and X. 1 ,, a function of X such as i// can be 
represented by a parabola, where 


= AX^ + BX + C 


(34) 


The parabola and corresponding points are shown in figure 6. The expressions for 
^i 1 j’ ^i j’ anc * *^i+l j can wr ^ en directly from equation (34) as 

tlJ = Ax2 j+ Bx u+ c 

<'l + l,j = AX ? + l,j +Bx i + l,j +c 
From figure 6 the following equations can be written: 


(35) 

(36) 

(37) 


X i,j~ X i-l,j + - 


x i + l,r x i-l,y+3 


Substituting these expressions into equations (36) and (37) gives 


^i,j = A ( X l-M + 5) +B ( X i-l.i + 5) + C 


(38) 


*i + l f j“ A ( X i- lf j + *J +B ( X i-l,J + ^ +C 


(39) 


Therefore, given equations (35), (38), and (39), it is possible to solve for A, B, and C 
to give 
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(40) 


\ 


A = 0. 5 d (\L/. + . - 2 il/. . + \b. - .) 


B = 0. 5 d(4^ j - 3l#Vlj , - ^. +1; .) - 2AX 1 _ 1> j 

C = " BX i-l, j " Axf-l, j 

The following derivatives can be written by using equations (40) to (42): 


(41) 

(42) 


#i_l 


— = 2 AX. 1 , + B = B' = 0. 5 d(4ii/. , - 3^. , . - \b. . .) 
0X. , . 1_1 ’J bJ ^ 1 - 1,1 v i+i,r 


i-1, j 


(43) 


d\f/ i . 

— it* = 2 AX, • + B = 2A 
9X i,J 


(X i, . + l\ + B = 2 AX. t +B+M =B ’+M 
\ 1 A,J d/ l_1 » J d d 


dxf/. ■ 


b3 


(44) 


= 2AX. , ; + B = 2A|X. 

® 1 + x,i 1+1,1 


( X i- 1 , j + f) + B = °- 5 *id,i - 4 *i, j + V'l-l.j) < 45 > 


...2 ...2 ...2 1+1 .J Vl - J 


ax? n . ax? . ax? , 

1-1,3 bJ i+i, J 


(46) 


Having developed these basic equations with the dummy variable i//, we can now 
proceed to develop the various terms in equation (33) by using the finite difference for- 
mat developed in this section. The following equation is written for the point (i, j): 


— /4 h 1 / 2 — ^ = 0. 5 d 

3X\t 7 ax 


>u 1+,J W i+1>j li-M l " ,J Wi-u 


(47) 


By using equations (42) and (44) the following equations can be written: 


3h\ 


ax) =°' 5 d( 3 H i + U- 4 H t,j + H i-M) 
' i+ 1, 3 


(48) 
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( 49 ) 


f) = °- 5 d( " H i+ij + 4H ij ~ 3H i-l, j } 
\ 3X/ i-l, j 


Substituting equations (48) and (49) into (47) gives 


—U H 1/2 ^ ) = 0. 25 d 2 Jh. - i (3H. . . - 4H. . + H. ^ .) 

ax\^ ax/ ^ i+1 j J 1+1,2 1+ ’•> 1,3 




(50) 


Thus, the following equation can be directly written: 


2-(jL = 0. 25 d 2 

ax\rj ax/ 


^i+k 


1 < s * i+ i, j - 4 * i, i + «i-i, j> + j - j + ®*i-i, j> 

^i+lj ^i-lJ 


(51) 


i-(pH) = 0.5d(p 1+ljj H i+1>r p 1 . 1;j H i . 1>j ) 


ax 


(52) 


The derivatives with respect to Y can be obtained directly by substituting c for d, 
the subscript i, j-1 for i-1, j, and the subscript i, j+1 for i+1, j. 


J_/£ H l/2iH\ o.25c 2 

SY\rj 3Y/ 


% j+1 


_a_/_£ a^\ _ o. 25 c 2 

ay ytj by) 


^VW H u + i- 4 H i,r 3 H u-!> 

j+1 C‘?^ - drh j. rh .) J. - bJ rl 


(53) 


(3« iJ+1 - 40^ + ^j.i) +z^ (0 iJ+1 - + 30- j.i) 

j+1 %j-l 


(54) 


^(pH) = 0.5 c(p l j+1 H ijJ+1 


(55) 
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ii i ill ■ ii 


Substituting equations (50) to (55) into equation (33) while collecting terms gives 


+ B i,Aj-i + c i,i*i-i,j + - hAj - M i,r° < 56 > 


where 




\r 3ll i* 1,1 + "n, j 
*i,j = (£f (( ‘i,i + i +3 "i,i-i ) 


C i, j "i+l, j + 3fi l- 1, j 

®i>j ’ (jj) 2 * 3 "i,j+l 

i i-i) 




(57) 

(58) 

(59) 

(60) 

(61) 


(3H i + l, J ‘ 4H U + W + M i-1» J V 5 ^ (H i + U - 4H U + + (i) 

J 

[“i.J+1 VVm (3H i,j»i - 4H U * H i,j-1> * "1,1-1 (h u*i - 4H i,i * 3H i,j-i)] 


\r 


24Ub 


^3/2 [ cos - ? i-i,j H i-i,j ) +£J ir j + i - "i,j-i H i,j-i ) ] 

J 

(63) 


Boundary Conditions 

The boundary conditions are the following: 

(1) At the edges of the rectangular computation zone (fig. 3) the pressure will be 

zero, therefore implying that <p is also zero. Specifically, this means that along the 

bottom of figure 3, 4> i l = 0; along the left side, 0. . = 0; along the top, 0. 20 = °; and 

along the right side, 0„,- . = 0. 

J 
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(2) At the cavitation boundary, P = 0 for 

3P = 3P _ o 
3X 3Y 

This condition is commonly known as the Reynolds condition and will be satisfied by 

simply resetting 0. . to zero whenever it occurs as a negative value. 

3 

(3) The term 0^ . is not allowed to exceed 4> max > where <P max = 60 H (where 0 H 
is equal to 0 when 0 is evaluated at the center of a Hertzian contact). 


Initial Conditions 

Outside the Hertzian contact ellipse the pressure is initially assumed to be zero, 
and therefore 0=0 there. That is, if P = 0, 


0 = 0 


when 


(X - m) 2 + (Y - l) 2 > 1 

Inside the Hertzian contact ellipse the pressure is initially assumed to be Hertzian; that 
is, 


or 


P = 


3F 
27rabE 


- Jl - (X - m) 2 - (Y - l) 2 


* - 3FhV2 Vl - (X - m) 2 - (Y - i) 2 


23:abE 


(64) 


(65) 


when 


(X - m) 2 + (Y - l) 2 < 1 
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Relaxation Method 


If subscript n is the iterant and . is the particular solution to be found, the 
relaxation method known as the Gauss-Seidel method can be expressed as 


^i,j,n+l ^i, j,n -X ^ 


- 3 " A i> n ' B i, A, j-l^n+1 ’ C i, j, n +1 ~ D i, j+1, n + 


( 66 ) 


L U 


0 . . 
h n 


where X is an overrelaxation factor which is initially set to 1. 4. Therefore, equa- 
tion (66) is used starting with node (2, 2) then continuing with (3, 2), . . . , until (M, 2); 
followed by (2, 3), (3,3), . . ., (M, 3); and ending with (2, N), (3, N), . . ., (M, N), 
where 


M = (n + m)d - 1 
N = 2lc - 1 


The relaxation procedure described by equation (66) is continued until 
N M 

I 

j=2, 3, . . . i=2 
where = 0. 1, initially. 


z 

,3,. 


K,j,n+1 ~ ^i,j,n 
^i, j,n+l 


<«i 


Dimensionless Pressure Loop 


The relaxation method provides values of j for every point within the mesh. 
Having determined 0. j, we can write the dimensionless pressure as 




(67) 


With these new values of the dimensionless pressure, new values of the dimensionless 
viscosity, ^density, and film thickness can be evaluated. Thus, the coefficients of equa- 
tion (66) (A^ j, B. j, . . . , j) should also change. Accordingly, it is necessary to 

reenter the relaxation loop. This process is continued until the following inequality is 
satisfied: 
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z 


z 


p - p 

i,j.n+l i,j,nl 
P i, 3, n+1 




NORMAL APPLIED LOAD 

The central dimensionless film thickness H q was initially estimated. The next 
task is then to find the correct value for H q . In order to do this, the integrated normal 
applied load must be evaluated, where 


/'m+n rll 

= E’ab J Q J 0 P dY dX 


( 68 ) 


By applying Simpson's rule, this double integral can be rewritten as 


v; 4E'ab 


2 (2 p 2, 2 

+ P 2, 3 

+ 2P 2,4 +P 2, 5 

+ . . 

• + 2P 2, 21-2 + P 2, 2 2-1 

( 2P 3,2 

+ P 3, 3 

+ 2P 3,4 +P 3, 5 

T . 

■ + 2P 3, 21-2 + P 3, 22 -1 

2 ( 2P m + n-2,2 

+ P m+n-2, 3 

+ 2P m+n-2, 4 + P m+n-2, 5 

+ . . 

' + 2P m+n-2, 21-2 + P m+n-2, 

^ 2P m+n- 1, 2 

+ P m+n-l, 3 

+ 2P m+n-l, 4 + P m+n-l, 5 

+ . . 

' + 2P m+n-l, 1-2 + P m+n-l, 




This equation can be rewritten simply as 

m+n-1 



i=2, 3, . . . 

where 

= 0 if i is odd 

q^ = 1 if i is even 
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q 3 = 0 if j is odd 

q 3 = 1 if j is even 

By using the trapezoidal rule, the normal applied load can also be expressed as 


m+n-1 

F - E ' ab 

cd / y 

1=2,3,. . 


2Z-1 




( 70 ) 


Both equations (69) and (70) will be programmed to serve as a check on one another. 

When the dimensionless film thickness is evaluated from equation (28), the smallest 

value of H, called H . , is retained as well as the location of H . . The location of 

min min 

H min is designated by X = X* and Y = Y*. The current values of the integrated nor- 
mal applied load F, the initially specified normal applied load F, the minimum film 
thickness H m ^ n , as well as the location of the minimum film thickness (X*, Y*), are 
used to make a new estimate for the central film thickness H q , which can be expressed 
as 



b 2 (X* - m) 2 _ a 2 (Y* - l) 2 _ w(X*, Y*) 

2RR 2RR R 

A y 


(71) 


where A is a constant, initially set at 0. 5. 

With the new value of the dimensionless central film thickness, the film thickness 
must be recalculated and reentry into the relaxation and pressure loop is required. 
This process is continued until the following inequality is satisfied: 



H 


o 


< ^3 


where is 0. 0005 initially. 


FLOW RATE 


The mass flow rate in the X and Y directions shown in figure 3 can be written as 
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X 12ij 8X 


<1^ 


= V ph -flS-ig, 


12)) 3Y 

By making use of equation (20), the preceding equations can be made dimensionless: 


Qx 


Q, 


^o 

_ / u v 

)pH 

pRH 3 

ap 

Po e ’ r2 

\REV 

12brj 

ax 

q Y r 'o 

II 

-3 

o 

I pH 

_ pRH 3 

ap 

Po E ’ r2 

\ RE V 

1 

12a?7 

3Y 

i 9 and 

U = Vrj 

o/ E 

’R, these ei 


as 




where 


(72) 


(73) 


(74) 

(75) 


ap 

ax 



p 


i-l, j 


) 


(76) 


ap 

3Y 




(77) 
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FLOW CHARTS 


Figures 7 and 8 are flow charts for the numerical solution on the digital computer of 
the equations developed in the analysis. Figure 7 is the flow chart of the main program. 
There are essentially three loops within the main program: In the relaxation loop, 

^i, j,n+l is generated. In the pressure loop, the new values of cj> i . n+1 of the relaxa- 
tion loop result in new values of pressure P^ j n+ ^ which in turn result in new values 
of film thickness j, new values of viscosity 77 . j, and new values of density .. 

The final loop is the central film thickness loop H q which ensures that the integrated 
normal applied load agrees with the initially specified value. 

Figure 8 shows the flow chart of the subroutine SUB 6 . Here it can be seen that a 
number of calculations occur only once and need not be repeated on reentry to this sub- 
routine. With a new pressure distribution the elastic deformation is recalculated and 
leaves this subroutine with a new film thickness and, therefore, a new d>. ■ . 


CONCLUDING REMARKS 

A procedure for the numerical solution of the complete, isothermal, elastohydro- 
dynamic lubrication problem for point contacts is outlined. This procedure calls for the 
simultaneous solution of the elasticity and Reynolds equations. In the elasticity analysis 
the contact zone was divided into equal rectangular areas. It was assumed that a uniform 
pressure was applied over each area. In the numerical analysis of the Reynolds equa- 
tion the parameter <p = PH^^, where P is dimensionless pressure and H is dimen- 
sionless film thickness, was introduced in order to help the relaxation process. 

The analysis is applicable to the complete range of elliptical parameters and is 
valid for any combination of rolling and sliding within a point contact. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, May 28, 1975, 

505-04. 
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APPENDIX - SYMBOLS 


A, B, C 

A,B,C,1 

D,L,MJ 

a 

a 

b 

b 

c 

D 

d 

E 

E' 


g 

F 

F 


G 

H 

H min 

H o 

h 

h min 

h o 

J 

k 

P 


constants defined in eq. (34) 
relaxation coefficients 

semimajor axis of contact ellipse 
a/2c 

semiminor axis of contact ellipse 
b/2d 

number of equal divisions of semimajor axis 
defined by eq. (30) 

number of equal divisions of semiminor axis 
modulus of elasticity 
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elliptical integral of second kind 

normal applied load 

integrated normal applied load 

elliptical integral of first kind 

dimensionless material parameter 

dimensionless film thickness, h/R 

dimensionless minimum film thickness, h min /R 

dimensionless central film thickness, h^R 

film thickness 

minimum film thickness 

central film thickness 

function of k (eq. (9)) 

elliptical eccentricity parameter, a/b 

dimensionless pressure, p/E T 
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p 

Q 

q 

R 

r 

U 

u 

V 
v 
W 
w 
w* 

x, X,xl 

y, Y, Y J 

z 

a,/3 

r 

v 

V 

% 

9 

\ 

P 

v 

P 

P 

<P 

<P 


pressure 

dimensionless mass flow rate, p?7 0 /p 0 E'R^ 
mass flow rate 
effective radius 
radius of curvature 

dimensionless speed parameter, V-q^RE' 
surface velocity in X direction 

Vu 2 ♦ v 2 

surface velocity in Y direction 

dimensionless load parameter, F/E'R R 

x y 

total elastic deformation 

total elastic deformation at center of contact 
coordinate systems defined in report 

viscosity pressure index, a dimensionless constant 

constants defining fluid used (eq. (23)) 

curvature difference 

lubricant viscosity 

dimensionless viscosity, r]/r\ 0 

atmospheric viscosity 

angle of resultant velocity vectors (eq. (18)) 
overrelaxation factor 
p / V 

Poisson ratio 
lubricant density 
dimensionless density, p/p Q 
defined by eq. (75) 

PH 2 / 2 
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Subscripts: 

A solid A 

B solid B 

x, y coordinate system defined in report 
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(a) Martin conditions: rigid solids and 
isoviscous lubricant. 


(b) Hertzian conditions: elastic solids and dry 


■M 


1— b —I 


tc) Elastohydrodynamic conditions-, elastic solids and 
Newtonian lubricant. 

u* . ^run . 1. 6 G ft 6 U°- 7 
R w 0.13 

Figure 1. - Lubrication and contact conditions and 
minimum-film-thickness formulas. 


Figure 2. -Geometry of contacting elastic solids. 
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Figure 3. - Division of area in and around contact zone into 
equal rectangular areas. 
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Figure 8. - Flow chart of subroutine SUB6. 
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Figure 7. - Flow chart of main program. 
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